#include "itgr_eq_sol.hpp"
#include <cmath>
#include <iostream>
using namespace std;

double k(double t, double s) { return std::exp(-s - t); }

double h(double x) { return std::exp(x) - std::exp(-x); }
/* f(x) - \int_0^1 k(x,y)f(y)dy = h(x) */
//精确解 exp(x)

double k2(double s, double t) { return 1 / (s + t + 1); }
double h2(double t) {
  return t - 1 + (t + 1) * std::log((2 + t) / (1 + t));
} // 精确解 x

//double k3(double t, double s) {
  //if (s > t)
    //return 2 * t * (1 - s);
  //else
    //return 2 * s * (1 - t);
//}
//double h3(double x) { return -2 * std::cosh(x); }

double k4(double x, double y) { return std::exp(-2 * x * y); }
double h4(double x) {
  return exp(-x) + exp(-2 * x - 1) / (2 * x + 1) - 1 / (2 * x + 1);
} // exp(-x)

int main(int argc, char *argv[]) {
  Gauss_IE_K_Matrix gkm(9);
  //auto k_m = gkm(k2, 0.0, 1.0, 3);

  //auto sol = gauss_ie_sol(k_m, h2);

  //for (int i = 0l; i < k_m.x_list.size(); ++i) {
    //std::cout << sol[i] << "    " << (k_m.x_list[i]) << std::endl;
  //}
  auto k_m = gkm(k4, 0.0, 1., 3);

  auto sol = gauss_ie_sol(k_m, h4);

  for (int i = 0l; i < k_m.x_list.size(); ++i) {
    std::cout << sol[i] << "    " << exp(-k_m.x_list[i]) << std::endl;
  }

  return 0;
}
